      program SCAR calculation ! Drew Shindell, Oct 23, 2013
      parameter (nf=8) ! number of forcers
      real*8 delT(0:1999),dT(nf,1999),iyr(0:1999),rf(nf,1999),
     & dam(nf,1999),disc(5,nf,1999),disc_sum(5,nf),emiss(0:1999),
     & dTref(0:1999),dTsum(nf,1999),rf_cc(nf,1999),dTc(nf,1999),
     & disc_tot(5,nf),dam_hc(5,nf),rng(5,nf)
      character*4 char4,fcr(nf)
      real*4 rs
      real*8 CH4(2,14999),rdisc(5),rcc_aer(5),dam_cl(5,nf)
      real*8 C,cd,rf_cc_CH4,
     * f,co2_in,ch4_in,n2o_in,co2_new,ch4_new,n2o_new,Fco2,Fch4,Fn2o
      real*8 CO2(0:14999),sco2(0:14999),e(1:nf)
      logical carbcyc
      data fcr/'CO2','CH4','BC','SO2','CO','OC','N2O','HFC'/
      data rdisc/1.05,1.03,1.014,1.025,1.04/
      data rcc_aer/0.91,0.79,0.61,.73,.73/

c     Use below if need to recalculate rcc coeffs (set to 1, then run w & w/o c-c)
c      data rcc_aer/1.,1.,1.,1.,1./

c     Native unit Mt emission
      un=1.E-6 !1.E-9 would be kg, 1.E-6 gives values per ton
      iyear =350 ! integration period, w/backstop dTlimit not important past ~150 yrs
      iyrplus=00 !gives SCC that many years in future rel to 2010
      irates=5 !number disc rates to output

c     Include c-c feedback for non-CO2 emissions
      carbcyc=.true.
c      carbcyc=.false.

c     **** Set up RF timeseries for each forcer (order in fcr array) ****
      rf(:,:)=0.

c     Carbon cycle section
C     integrate sources over time to calculate CO2 (source in GtC)
        CO2 = 0.
        co2rfsum=0.
        do j=1,iyear

c        Use the next two lines to test RF from pulse vs time horizon
         sco2(j)=un*1.E-3*12./44. !1 ton CO2 as conversion to ppm is for GtC
         if(j.ge.2)sco2(j)=0.

         do k=1,j
          cd = cdecay(real(j-k))
          CO2(j)=CO2(j)+sco2(k)*cd
         enddo
         CO2(j)=CO2(j)*0.47 ! roughly 0.47 ppmv/GtC (physical constant)
c        set base CO2
         CO2_init=260.d0 ! 2000 CO2 (ppmv)

         CO2(j) = CO2_init + CO2(j)
         rf(1,j)=5.3d0*log(co2(j)/CO2_init)
         if(j.le.100)co2rfsum=co2rfsum+rf(1,j)
        end do

C     Methane
      tau=12.4 ! Perturbation lifetime
      ch4rfsum=0.
      do j=1,iyear
        rf(2,j)=un*.00032*(exp(-real(j-1)/tau))
        if(j.le.100)ch4rfsum=ch4rfsum+rf(2,j)
      enddo
c     compare integrated 100 yr RF for CO2 and CH4
      print*,"GWP100 CH4 = ",ch4rfsum/co2rfsum

c     Set BC RF (1 year only, all others zero)
      rf(3,1)=un*(-0.40/(1.184-5.339)) ! based on UNEP modeling, W/m2 per Tg BC
c     UNEP had BC change for all measures case of 5.339 Tg -> 1.184 Tg, RF was -0.40

c     Set SO2 RF (1 year only, all others zero as virtually all RF via SO4 and NO3)
c     Shindell et al., Science, 2009 had 131 Tg SO2 and RF = -0.869 incl AIE
      rf(4,1)=un*(-0.869/131.0)

c     Set CO RF (1 year part via ozone and aerosols, CH4 tau part as well)
c     Shindell et al., Science, 2009 had 662 Tg CO and RF = 0.224 from O3, SO4, NO3 (incl AIE)
c      also was 0.080 RF from CH4
      rf(5,1)=un*0.224/662.0
      do j=1,iyear
       rf(5,j)=rf(5,j)+un*(.080/662.)*(exp(-real(j-1)/tau))
      enddo

c     Set OC RF (1 year only, all others zero)
      rf(6,1)=un*0.15/(2.632-13.933) ! initial estimate based on UNEP modeling, W/m2 per Tg OC
c     UNEP had OC change for all measures case of 13.933 Tg -> 2.623 Tg, RF was 0.15

c     N2O (based on its GWP and lifetime)
      rfsumn2o=0.
      do j=1,iyear
        rf(7,j)=un*(exp(-real(j-1)/121.0))*0.000488
         if(j.le.100)rfsumn2o=rfsumn2o+rf(7,j)
      enddo
c     compare integrated 100 yr RF for CO2 and N2O
      print*,"GWP100 N2O = ",rfsumn2o/co2rfsum

c     HFC-134a (based on its GWP and lifetime)
      rfsumhfc=0.
      do j=1,iyear
        rf(8,j)=un*(exp(-real(j-1)/13.4))*0.0118 
         if(j.le.100)rfsumhfc=rfsumhfc+rf(8,j)
      enddo
c     compare integrated 100 yr RF for CO2 and HFC-134a
      print*,"GWP100 HFC-134a = ",rfsumhfc/co2rfsum

C     ***** Global temperature response *****

       rs=(3.75/3.75)*(1./1.07) !clim sens * (1/Had CM clim sens)

C     GTP impulse response fcn from Boucher et al., ERL, 2009
c       calculate impulse-response fcn
        do j=1,iyear
          delT(j)=(rs*0.631/8.4*EXP(-(j)/8.4)+
     &     rs*0.429/409.5*EXP(-(j)/409.5))
        enddo

c     define reference temperature trend
      dTlimit=4.47
      dTref(1)=0.8
      dTref(1)=dTref(1)+(0.6/40.)*iyrplus !gives SCC that many years in future
      dTref(0)=0.
      dTrate=.015 ! .6 C is roughly 1970-2010 observations, .6/40=.015
      do j=1,iyear
c       Include increase at first, then slowdown past threshold assuming societal actions kick in
        dTref(j+1)=dTrate*(1.013**j)+dTref(j)
        if(dTref(j-1).ge.4)then
         dTrate=.008
         dTref(j+1)=dTrate+dTref(j)
c        max dT limit
         if(dTref(j+1).ge.dTlimit)dTref(j+1)=dTlimit
        endif
      enddo

      dT(:,:)=0.
c     Calculate global temperature timeseries
        do j=2,iyear ! j loops through all years
         do k=1,j-1 !k is starting year through previous year of j loop
           do i=1,nf ! forcers
            dT(i,j)=dT(i,j)+delT(j-k)*rf(i,k)
           enddo
         enddo
        enddo

c      include c-c feedback for non-CO2 forcers
       if(carbcyc)then
        rf_cc(:,:)=0.
        rf_cc_CH4=0.
c       calculate forcing due to c-c feedback
        do j=2,iyear ! j loops through all years
         do k=1,j-1 !k is starting year through previous year of j loop
          do i=2,nf    
c          rf(1,j) is fcg at yr j due to 1 Mt CO2 in yr 1
c          extra CO2 emitted is 1 GtC per K (44/12*1.E12 in kg CO2); Collins et al., ACP, 2013
           rf_cc(i,j)=rf_cc(i,j)+(44./12.)*1.E12*dT(i,k)*rf(1,j-k)*
     &8.E-14/co2rfsum
          enddo
         enddo
         if(j.le.100)rf_cc_CH4=rf_cc_CH4+rf_cc(2,j)
        enddo
c       add effect of c-c feedback forcing
        dTc(:,:)=0.
        do j=2,iyear ! j loops through all years
         do k=1,j-1 !k is starting year through previous year of j loop
           do i=2,nf ! forcers
            dTc(i,j)=dTc(i,j)+delT(j-k)*rf_cc(i,k)
           enddo
         enddo
          do i=2,nf
           dT(i,j)=dT(i,j)+dTc(i,j)
          enddo
        enddo
c      compare integrated 100 yr RF with c-c fdbk CH4
       print*,"GWP100 CH4 w cc-fdbk = ",(ch4rfsum+rf_cc_CH4)/co2rfsum
       else
       endif

      dTsum(:,1)=0.
C     add reference and dT
      do j=1,iyear
        do i=1,nf ! forcers
          dTsum(i,j+1)=dT(i,j)+dTsum(i,j)
          if(dTref(j).ge.dTlimit)then
            dT(i,j)=dTlimit
          else
            dT(i,j)=dT(i,j)+dTref(j)
          endif
        enddo
      enddo

c     ***** Calculated Damages *****

c     Calculate global discounted damage function (% of GDP) due to climate (dT)
      disc_sum(:,:)=0.
      disc_tot(:,:)=0.
        do j=2,iyear ! j loops through all years
c        5th discount rate is declining with time (DDR)
c        use 250 yr e-folding approx fm Arrow et al, Science, 2013.
         rdisc(5)=1.+0.01*(4.*exp(-real(j)/250.))

          do i=1,nf ! forcers
c     Damages based on .00228=1.8%GDP/(2.5C)^2 econ damage coefficient
c     Assume GDP grows 2% per yr increase - this gives 2100 GDP of 355 trillion, as in IWG
          dam(i,j)=.00288*(dT(i,j)**2-dTref(j)**2)*(1.02**(j+iyrplus-1))
c     test function in Ackerman & Stanton fm DICE, nearly same results
c       dam(i,j)=(1./(1.+(dT(i,j)/18.8)**2))-(1./(1.+(dTref(j)/18.8)**2))
c       dam(i,j)=dam(i,j)*(1.018**(j-1)) 
          do id=1,irates
           disc(id,i,j)=dam(i,j)/(rdisc(id)**(j-1))
          enddo

           disc_sum(:,i)=disc_sum(:,i)+disc(:,i,j)
          enddo
        enddo

c       Convert from damages in %GDP to SCAR in $2007
c       use actual worldwide GDP to get dollars (50 trillion 2010 GDP from IWG)
        rfactor=50.E12

        do i=1,nf ! forcers
          disc_sum(:,i)=disc_sum(:,i)*rfactor
c         valuation for non-CO2 enhanced vs CO2 as others more damaging to plants since 
c         no offset by CO2 fertilization; use 10% (range 5-20%)
          if(i.ne.1)disc_sum(:,i)=disc_sum(:,i)*1.1
        enddo

        write(*,'(15x,a5,8(5x,a4))') 'rate ',fcr(:)
        do i=1,irates
         write(*,'(a8,1x,f10.3,8(1x,i8))') 
     &    "Climate ",rdisc(i),int(disc_sum(i,:)+0.5)
         disc_tot(i,:)=disc_tot(i,:)+disc_sum(i,:)
         dam_cl(i,:)=disc_sum(i,:)
        enddo

        write(*,*)
c       aerosol climate valuation due to regional change

c       Don't use full regional aerosol val as part of impact from C-C response rather than aerosols
c       run twice, calculate fraction due to c-c response, remove from regional (rcc coeffs)

        do i=1,irates
         write(*,'(a8,1x,f10.3,8(1x,i8))') 
     &    "Reg Aer ",rdisc(i),0,0,
     &    int(disc_sum(i,3)*1.6*rcc_aer(i)+.5),
     &    int(-disc_sum(i,4)*(1.6*rcc_aer(i)+1.)+.5),0,
     &    int(-disc_sum(i,6)*(1.6*rcc_aer(i)+1.)+.5),0,0
        enddo
        
        disc_tot(:,3)=disc_tot(:,3)*(1.+1.6*rcc_aer(:))   !BC, all postive=2.6 total
        disc_tot(:,4)=-disc_tot(:,4)*(1.6*rcc_aer(:))  !SO2, 0.5-0.5*4.2=-1.6 
        disc_tot(:,6)=-disc_tot(:,6)*(1.6*rcc_aer(:))  !OC as SO2
        disc_sum(:,:)=0.

c       Estimate health costs of climate change based on WHO estimate of additional premature deaths
c        due to climate change (current) of 0.8 C gives 150,000
        do j=2,iyear ! j loops through all years
         rdisc(5)=1.+0.01*(4.*exp(-real(j)/250.))
         do i=1,nf 
c       varies with dT^2 with 150000/(0.8 degree) at $7.4 million/life * ratio country-specific VSl/uniform 
c         VSL from UNEP; divide by 1E6 to get per ton rather than per Mt
          dam(i,j)=(dT(i,j)**2-dTref(j)**2)*
     &     ((150000./0.8)*7.4/un*3.46/14.8)*0.995**(j+iyrplus-1)
c          population incr w time at 0.4% (gives 9 billion in 2100), but baseline mortality decreases
c           here assume net is .5% decrease

          do id=1,irates
           disc(id,i,j)=dam(i,j)/(rdisc(id)**(j-1))
          enddo
c          Climate-health based on net (clim+reg aer) so * (clim+reg aer)/clim
           if(i.eq.3.or.i.eq.4.or.i.eq.6)
     &      disc(:,i,j)=disc(:,i,j)*disc_tot(:,i)/dam_cl(:,i)
          disc_sum(:,i)=disc_sum(:,i)+disc(:,i,j)
         enddo
        enddo

        write(*,*)
        write(*,'(15x,a5,8(5x,a4))') 'rate ',fcr(:)
        do i=1,irates
         write(*,'(a8,1x,f10.3,8(1x,i8))') 
     &    "HlthCli ",rdisc(i),int(disc_sum(i,:)+.5)
         disc_tot(i,:)=disc_tot(i,:)+disc_sum(i,:)
        enddo
        dam_hc(:,:)=disc_sum(:,:) !save health-climate damages for later use in totals
        disc_sum(:,:)=0.

C       Calculate damages due to health impacts of direct composition changes
c       Methane impacts are 1080 +- 721 per tonne from Shindell et al., Sci, 2012.
c         This value for methane is using country-specific VSL.
c       This is total integrated future value with no discounting (as it's equil response)
c        so normalize by 1080/11349 (denom is sum of 1080 with decay)
        do j=1,iyear
         rdisc(5)=1.+0.01*(4.*exp(-real(j)/250.))
c        Health impact proportional to ozone response, which follows CH4 RF in time
c        also account for 2010 population is 6.5 billion, 2030 is 8.4 billion
c        and for VSL, 7.4 million $ is 2006, 9.5 used in 2030
c        assumes no change in 2030 vs 2010 baseline mortality (as in original calc
c        described in Shindell et al/UNEP as no mort proj available)
         rfactor=(rf(2,j)/rf(2,1))*(7.4/9.5)*(6.5/8.4)*1080./11349.
         do id=1,irates
           disc(id,2,j)=1080.*rfactor/(rdisc(id)**(j))
         enddo
         disc_sum(:,2)=disc_sum(:,2)+disc(:,2,j)
        enddo

c       calculate aerosol health impacts
c       UNEP/Science paper used linear PM2.5 health impact, as did NCC 2011 paper
c       so take fraction per aerosol type * total PM2.5 impact and
c       divide by total emissions to get impact per unit emission
c       total current PM2.5 impact (NCC paper): 2,000,000 lives
c       valuation is UNEP 2030 country-specific, * ratio 2030/2000 pop, *ratio UNEP/PD deaths =$3.46trillion
c       scale to match new 2013GBD est of 3.2 million lives ($3.46->$5.54 trillion)

c       fractional contrib to PM2.5: SO4 = 37%, BC = 5.5%, OC = 32% (NO3 19%, SOA 6%) (w pop weighting)
c       total emissions: 121.8 Tg/yr SO2, 8.9 Tg/yr BC, 64.7 OC (anthro)
c       use valuation over 1E6 to convert from Tg to tons reduced
        disc(:,3,1)=5540000.*0.055/8.9 ! BC
        disc(:,4,1)=5540000.*0.37/121.8 ! SO2
        disc(:,6,1)=5540000.*0.32/64.7 ! OC
        disc_sum(:,3)=disc_sum(:,3)+disc(:,3,1)
        disc_sum(:,4)=disc_sum(:,4)+disc(:,4,1)
        disc_sum(:,6)=disc_sum(:,6)+disc(:,6,1)

        write(*,*)
        do i=1,irates
         write(*,'(a8,1x,f10.3,8(1x,i8))') 
     &    "HlthCmp ",rdisc(i),0,int(disc_sum(i,2)+.5),
     &    int(disc_sum(i,3)+.5),int(disc_sum(i,4)+.5),0,
     &    int(disc_sum(i,6)+.5),0,0
c        Uncertainty range due to 50% health-climate impact plus 80% on health-comp impact
c        80% fm # premature mortalities/health damages in UNEP
         rng(i,:)=sqrt((disc_tot(i,:)*.5)**2+(disc_sum(i,:)*.8)**2)
         disc_tot(i,:)=disc_tot(i,:)+disc_sum(i,:)
        enddo
        disc_sum(:,:)=0.

C       Calculate damages due to agricultural impacts of direct composition changes
c       Methane impacts are 29 +- 8 per tonne from Shindell et al., Sci, 2012.
        do j=1,iyear
         rdisc(5)=1.+0.01*(4.*exp(-real(j)/250.))
c        Crop impact proportional to ozone response, which follows CH4 RF in time
         rfactor=(rf(2,j)/rf(2,1))*1080./11349.
         do id=1,irates
           disc(id,2,j)=29.*rfactor/(rdisc(id)**(j))
         enddo
         disc_sum(:,2)=disc_sum(:,2)+disc(:,2,j)
        enddo

        write(*,*)
        print*,"CH4 ag due to comp at 5%",
     &   int(disc_sum(1,2)+.5)
        print*,"CH4 ag due to comp  1.4%",
     &   int(disc_sum(3,2)+.5)
         disc_tot(i,2)=disc_tot(i,2)+disc_sum(i,2)

c       write totals
        write(*,*)
        write(*,'(15x,a5,8(5x,a4))') 'rate ',fcr(:)
        do i=1,irates
         write(*,'(a8,1x,f10.3,8(1x,i8))') 
     &    "Tot kg ",rdisc(i),int(disc_tot(i,:)-dam_hc(i,:)*0.5+.5)
        enddo
        write(*,*)
        do i=1,irates
         write(*,'(a8,1x,f10.3,8(1x,i8))') 
     &    "Range  ",rdisc(i),int(rng(i,:)+.5)
        enddo

      end

c     Functions

      function cdecay(t)
      real*4 t ! time in years
      parameter (n=3)
      dimension a(0:n),tau(0:n)
c**** Exponential values taken from Bern Carbon Cycle Model
c****  (Siegenthaler and Joos, 1992)
c     IPCC AR4 c-c version (Forster et al 07) as quoted in Collins et al., ACP, 2013
      data a/0.217,0.259,0.338,0.186/
      data tau/100000.,172.9,18.51,1.186/

      cdecay = 0.
      do i=0,n
        cdecay=cdecay+a(i)*exp(-t/tau(i))
      end do

      return
      end

      real*8 function f(M,N)
      real*8 M,N

      f=0.47*log(1.+2.01d-5*(M*N)**(0.75d0)+5.31d-15*M*(M*N)**(1.52d0))
      return
      end

